function [Pxx,f] = myWelch(x, fs, windowLength, overlap)
% x: 输入的信号序列
% fs: 采样频率（单位：Hz）
% windowLength: 窗函数的长度（点数）
% overlap: 相邻段之间的重叠点数

% 计算每段的不重叠部分长度
noverlap = windowLength - overlap;
% 计算分段数量
nfft = windowLength;
% nfft = 2^nextpow2(windowLength); % 为了提高计算效率，取大于等于windowLength的2的幂次
numSegments = floor((length(x) - overlap) / noverlap);

% 初始化功率谱估计结果矩阵
Pxx = zeros(nfft, 1);

% 对每一段信号进行处理
for i = 1:numSegments
    startIndex = (i - 1) * noverlap + 1;
    endIndex = startIndex + windowLength - 1;
    segment = x(startIndex:endIndex); % 取出一段信号
    window = hanning(length(segment));
    pow_win = mean(abs(window).^2);
    windowedSegment = segment.* window; % 加汉宁窗（可更换其他窗函数）
    segmentFFT = abs(fft(windowedSegment, nfft)).^2/pow_win; % 计算该段的FFT幅值平方
    % segmentFFT = abs(fft(windowedSegment, nfft)).^2; % 计算该段的FFT幅值平方, 不对窗函数功率归一(U)
    Pxx = Pxx + segmentFFT; % 累加功率谱
end

% 平均功率谱
Pxx = Pxx / numSegments;

% 生成对应的频率轴
f = (0:(nfft - 1)) * fs / nfft;
end